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Abstract — In this article, we develop a continuous implicit seventh-eight method using interpolation and 
collocation of the approximate solution for the solution ofy’ — f(x, y) with a constant step-size. The method 
uses power series as the approximate solution in the derivation of the method. The independent solution was 
then derived by adopting block integrator. The properties of the method was investigated and found to be 
zero stable, consistent and convergent. The integrator was tested on numerical examples ranging from linear 
problem, Prothero-Robinson Oscillatory, Growth Model and Sir Model. The results show that the computed 
solution is closer to the exact solution and also, the absolutes errors perform better than the existing methods. 
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I. INTRODUCTION 

Numerical analysis predates the invention of modern 
computers by many centuries. Linear interpolation was 
already in use more than 2000 years ago. Many great 
mathematicians of the past were preoccupied by numerical 
analysis, and this is obvious from the names of important 
algorithms like Newton's method, Lagrange interpolation 
polynomial, Gaussian elimination and Euler's method. 
Despite the long and rich history of numerical analysis, 
"modern" numerical analysis is characterized by the 
synergy of the programmable electronic computer, 
mathematical analysis, and the opportunity and need to 
solve large and complex problems in applications. The need 
to advance in applications such as ballistics prediction, 
neutron transport, and non-steady multidimensional fluid 
dynamics motivated the development of the computer 
which depended strongly on advances in numerical analysis 
and mathematical modelling. 

1.1 Aim and Objectives 


The aim of this research work is to develop continuous 
method via fractional-step methods by using power series 
polynomial as basic functions for solving IVPs in ODEs. 
The objectives are to. 

i. develop hybrid block methods through continuous 
implicit seventh-eight method 

ii. analyzing the basic properties of the method; zero 
stability, order, convergence and region of absolute 
stability 

iii. implementing a self-staring implicit methods 

without the rigor of developing predictors 

iv. developing computer codes to implement the 

method derived using MATLAB 

v. examine the exact solution and the computed 

solution, if they are in agreement 

vi. compute the time taken for computation in seconds 

1.2 Significance of the Study 

The study is significant for the following reasons 

1. This work centered on the numerical solution of 
first order ordinary differential equations of initial 


http://www.aipublications.com/iicmp/ 


Page | 39 






Open Access 


International journal of Chemistry, Mathematics and Physics (IJCMP) 

[Vol-4, Issue-3, May-Jun, 2020] 
h ttps://dx. doi. org/10.22161/ijcmp. 4.3.2 

ISSN: 2456-866X 


value problem. Hence, the differential equations 
whose order is more than one are not considered. 

2. The proposed method makes use of power series 
polynomial, other polynomials are not considered 
in this work. 

3. Throughput this research work, linear fractional- 
step method is considered. Non hybrid points are 
not considered. 

4. Most of the modeled problems in ordinary 
differential equations do not have analytical or 
theoretical solution, hence the use of numerical 
solution is very important 

II. LITERATURE REVIEW 

This work is based on numerical solution of first order 
initial value problems of ordinary differential equation, 
which is of the form 

y'= f(x,y), y(a) = rj, a<x<b (1) 

where f is continuous with the interval of integration. It is 
assumed that f satisfies Litchitz condition which guarantees 
the existence and uniqueness of solution (1). 

Countless real life problems in sciences, engineering, 
biology and social sciences are model of first order ordinary 
differential equations. For instance. Environmental Models 
such as SIR Model, Growth model. Mixture model just to 
mention few. Olver (2008) and Kandasamy el al. (2005), 
reported that equation (1) is used in simulating the growth 
of populations, trajectory of a particle, simple harmonic 
motion, deflection of a beam etc. With the advent of the 
modern high speed electronic digital computers, the 
numerical integrators have been successfully applied to 
study problems in mathematics, engineering, computer 
science and physical sciences such as physics, biophysics, 
atmospheric sciences and geosciences, Jain et al. (2007) 

Many numerical integration schemes to generate the 
numerical solution to problems of the form (1) have been 
proposed by several authors. The one-step methods include 
those developed by Fatunla (1994), Aashikpelokhai (1991), 
Ayinde and Ibijola (2015), and Kama and Ibijola (2000), 
Ibijola (1998), Ogunrinde (2010) etc . These methods were 
constructed by representing the theoretical solution y(x) to 
(1) in the interval, \x n , X n+1 ] n> 0 by linear and non-linear 
polynomial interpolating functions. 


Furthermore, authors like Butcher (2003), Sunday et al. 
(2015), Lukman and Olaoluwa, Areo et al. (2011), Ibijola, 
Skwame and Kumleng (2011), Olangegan et al, (2015), 
Areo and Omojola (2017), Adeniran and Longe (2019), 
Adeniyi et al. (2006), Onumanyi et al, (1994), Sunday and 
Odekunle (2012), Badmus (2014), Jator (2014), Kayode and 
Adeyeye (2013) and Ibijola et al. (201) have all proposed 
linear multistep methods (LMMs) to generate numerical 
solution to (1). These authors proposed methods in which 
the approximate solution ranges from power series, 
Chebychev’s, Lagrange’s, Lucas polynomial, Legendre 
Polynomial and Laguerre’s polynomials. In spite of the 
improvement made by diverse researchers in developing 
schemes that solve (1), the order of the schemes are low 
resulting in lower accuracy when used to solve problems 
which is a setback. Hence, the setbacks by previous works 
motivated a need to develop a seventh-eight method for 
solving first order ordinary differential equations unlike 
other methods which is one-step or multi-step. The method 
will be analysed and test on some numerical examples to 
test the accuracy and efficiency of the new method in order 
to ascertain how reliable it is. 

III. MATERIALS AND METHODS 

The power series of the form (2) is considered 

S+(7~ 1 

y(*)= Z 

n =0 

( 2 ) 

The first derivative of (2) is given by (3) 

£ +cr- l 

/(*)= Z a nP'M )’ 

n =1 

(3) 

as an approximate solution to equation (1) where 
p n (x) = x n ,. 8 and cr are the number of distinct 

collocation and interpolation points respectively. 
Substituting the first derivative of (3) into (1) gives 

c+i - 1 

f(X,y(x),)= ^ a nP'n( X )’ 

n =1 

(4) 
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Collocating (4) at X — 0 




\° j 


— and interpolating (3) at X — — leads to a system of nine equations which is solved with the aid 

8 4 


of Maple Mathematical Software version 18.0 to obtain Cl n , n = 0 , 1 , •••,8 . The <2 n ' S obtained are then substituted into (3) to 
obtain the continuous form of the method 


y(x) = a 3 (t)y 3 +h 

— n -\— 

4 4 

where CC n and are continuous coefficients. 


A>( 0 /„+A( 0 / i +M)f !+y? 3 ( 0 / 3 + M)f 1 +A( 0 / 


n -\— — n -\— 

8 4 4 


n -\— — iu — 

8 2 2 


+A(0/ 3 +A( r )/ 7 

— n -\— — n -\— 

4 4 8 8 


(5) 


a 3 (f) = l 

T 


[3 0 (0 = —30240” ( 6144/?4 ^ “ 12288 h 3 t 3 + 8192 h 2 t 2 — 2048 h t 
+ 123) {4ht- l) 2 (4 /? / — 3 ) 2 

J3, (0 =^8o" ( 86016 h * f ~ 2 °2 7 52 h 5 t 5 + 176128 h A t A - 69312 h 3 ? 

¥ 

+ 11328 h 2 t 2 — 216 /? ? — 81) (4/z/-3) 2 

J3 j (0 = —(229376 h 6 t 6 — 507904 h 5 t 5 + 399360 h A t A - 133120 h 3 t 3 
¥ 

+ 15728/t 2 / 2 + 24 ht + 9) [4ht-3] 2 

P At) = -j— (215040 h 6 t 6 — 445440 h 5 t 5 + 317440 h A t A - 92448 h 3 f* 

_ 1890 

8 

+ 9528/?V — 136/tf —5l) (4 /? / — 3) 2 

J3 j (0 = - 3Q j, 4Q (3440640 h 6 t 6 — 6635520 h 5 t 5 +4311040 h 4 t A 
~2 

- 1155072 h 3 ? + 118032 /t 2 t 2 + 216 /z t + 81) (4 /z t - 3 ) 2 

/3 5 (0 = (28672 h 6 t 6 - 51200 h 5 t 5 + 30720 h 4 t A - 8000 /? 3 ? 

¥ 

+ 736 h 2 t 2 — 24 h t — 9) {4ht~3) 2 
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p 3 (/) = - 3Q ^ 4Q {4ht-3) (2752512 h 1 1 1 — 6586368 /z 6 t 6 


+ 5955584 h 5 t 5 - 2586624 h 4 t 4 + 545088 h 3 t 3 - 47696 h 2 t 2 - 492 h t 


- 369) 

/3 ? (/)=^V (768 /z 4 / 4 — 1152 /? 3 + 640 A 2 t 2 — 156 h t + 15) (4 ^ 

¥ 

— 3 ) 2 

The continuous method (5) is used to generate the main method. That is, we evaluate at X = X 7 

n+- 


y 


z y 3 + 


1 


n + L - y n + A. 967680 

8 4 


1375/^ + 123133/ j +41499/ , 

« + — n + 


11351/ , + 139849/ , - 88547/ 


2 4 

- 121797/ s 


n + ■ 


n + ■ 


n + ■ 


n + 


8 


+ 36799/ ' 


n + l 


(6) 


Continuous method (5) is also used to generate the additional method at the non-interpolation points. That is, we evaluate at 


X X 1 ’ X 1 ’ X 3' X 1 anc ^ X — X 


n -\— n -\— n -\— n -\— 

8 4 8 2 


5 

n -\— 


} ’ n y n + l~ 1120 


41/ w +27/ j +27/ , +216/ t +41/ 3 


n + ■ 


n + ■ 


n + • 


n + ■ 


V 


+ 272/ +216/ 

n + — 




(7) 


y 1 = y 3 + 


„ + _ „ + + 193536 

8 4 


55 / -4475/ . - 5805/ 


n + j 


n + J 


- 1871/ , - 1871 / 3 -4475/ 3 - 5805/ 5 + 55/ ? 

« + — w + — n + — « + — « + — 

4 


\ 
8 ", 


( 8 ) 


^ 1 => ? . 3 


« + — « + /- 3780 

4 4 


332/ j + 171/ , -4/ j +171/ 3 

W + T " + T w + ¥ M + T 


+ 612/ 3 +612/ s -4/ 


\ 


n + • 


n + • 


n + • 


(9) 
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y +3 =y 3 

n + — n + — 
° 4 


+ 


35840 


h 


n + ■ 


n + 


- 5535/ 5 +45/ / 

n + — n + — 


- 1685/ 3 - 2777/ 3 

« + T « + - 

5^-1153/ j +135/ j -40/ 

n + — n + — n - 

2 4 


7 1 3 

n + — n + — 
2 4 


+ 


1 


30240 




Z H V 

- 1363/ 3 -208/ 3 -4968/ s + 32/ 7 

w + — ,2 + — w + « + — 

o 




J 5 3 

n + — n + — 

4 


+ 


967680 


h 


63514 + 44797/ , + 11547/ 

« + — n 


o 4 y 

-2999/ , -47799/ 3 -26883/ 3 


« + ' 


101349/ 5 

n + — 


+ 1375/ ? 

« + — 


(ll) 


( 12 ) 


The methods derived in equation (6) to (12) are combined and implemented as a block in solving numerical problems of (1). 


IY. ANALYSIS OF THE METHOD 


4.1 Order and Error Constant of the Block Method 

We adopt the method proposed by Fatunla (1994), Lambert (1973) and Henrici (1962) to obtain the order of the method obtained 
from (6-12) as (8,8,8,8,8,8,8) T and Error constants as 


[-6.97e“ n , -5.47e“ u , -6.14e“ u , -5.62e -11 , -6.29e” n , -4.79e” u , -1.18e“ u ] r 

4.2 Zero Stability of the Block Method 

A block method is said to be zero stable, if the roots 


(13) 


det[/L4 (0, -A (,) ] = 0 


(14) 


of the first characteristic polynomial satisfy |/l| < 1 and for the roots with |/l| < 1, the multiplicity must not exceed the order of 
the differential equation. 
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f 1 

0 

0 

0 

0 

0 

0^ 


fo 

0 

0 

0 

0 

0 

0 



0 

1 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 



0 

0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 


z 

0 

0 

0 

1 

0 

0 

0 

- 

0 

0 

0 

0 

0 

0 

1 



0 

0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

0 

0 

1 



0 

0 

0 

0 

0 

1 

0 


0 

0 

0 

0 

0 

0 

1 



10 

0 

0 

0 

0 

0 

Ij 


10 

0 

0 

0 

0 

0 

1J 



A = z 6 (z-l) = 0, z = 0,0,0,0,0,0,0 

Hence the block is zero stable. 

4.3 Consistency of the Method 

Our new developed block method is consistent since the order of each of the method is greater than 1. 

4.4 Convergence 

According to Lambert, 1973, Brown 1977 and Chollom et al. (2007) the necessary and sufficient condition for a linear multistep 
method to be convergent is for it to be consistent and zero stable. 

Hence Convergence = Consistency + Zero-stability. 

Hence the new block method is convergent 

4.5 Region of Absolute Stability of the Proposed Method 

According to Ibijola et al. (2011), Areo and Omojola (2017), the stability matrix is expressed as 


M{z) = V + zB(M -zA)~ l U 

together with the stability function 

p(rj, z ) = det(r/I -M(z)) 

We express the block method obtained from (6-12) in the form 


Y 


A 

u 

hf(y ) 

Y i + 1 


B 

V _ 

Y i-l 


(16) 


(17) 


(18) 


where 
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1 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

: V-= 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 


0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

0 

0 

1 


139849 

4511 

123133 

88547 

1537 

11351 

275 

967680 

35840 

967680 

967680 

35840 

967680 

193536 

733 

71 

17 

1927 

13 

29 

1 

3780 

3360 

210 

30240 

420 

3360 

945 

1359 

1377 

5927 

3033 

1377 

373 

9 

7168 

35840 

35840 

35840 

35840 

35840 

7168 

181 

1 

223 

53 

1 

8 

1 

945 

35 

945 

3780 

35 

945 

945 

36725 

775 

4625 

13625 

1895 

275 

275 

193536 

21504 

21504 

193536 

21504 

21504 

193536 

27 

27 

17 

27 

27 

41 


140 

1120 

70 

1120 

140 

1120 


25039 

343 

20923 

20923 

343 

25039 

5257 

138240 

5120 

138240 

138240 

5120 

138240 

138240 


B := 


0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 


5257 


y n 


fn 



138240 

41 


y i 

n-\ — 

8 


f 1 

n-\ — 

8 



1120 

265 


y i 

n-\ — 

4 


/ i 

n -\— 

4 


y i 

n -\— 

8 

y 7 

n -\— 

8 

7168 

139 

Y = 

y 3 

n -\— 

8 

> f(y) = 

/ 3 
n -\— 

8 

. ^-i = 

y i 

n -\— 

8 


3780 


y i 


/ t 


[y n _ 


265 


2 


2 



7168 


y |5 


/,5 



41 




” 8 



1120 


y n i 


/, 7 



5257 


8 _ 


8 _ 




138240 


(19) 


The elements of the matrices A, B, U and V are substituted and computing the stability function with Maple software yield the 
stability polynomial of the method which is then plotted in MATLAB environment with the Newton-Raphson Method where 
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N = 400 and tol= 10 12 to produce the required absolute stability region of the method as shown in the figure 1 below. The 
graph is shows the region in which the method is absolutely stable. We conclude that the method is P-stable since the interval of 
periodicity lies between (0, 2.0). 


Region of Absolute Stability of the Block Method 



V. NUMERICAL EXPERIMENT 

This section deals with the implementation of the block 
method in solving initial value problems (IVP) of first order 
ordinary differential equations. The method is coded in 
MATLAB (R2012a) version environment using window 10 
as an operating system on Acer Laptop. The new developed 
method is tested on some problems to determine the 
performance of the new proposed methods. The error is 

defined as Error = |_y(x) — y n (x)| , where y(.l') is the 
exact solution, y n (x) is the computed result. 

5.1 Numerical Examples 

The following problems are chosen as test problems 
Problem 1: Consider the first order problem 


y' = —y, t(0) = l, 0<*<1, h = 0.1, 


( 20 ) 


Exact solution: y(x) = e 

Problem 2: Consider the Prothero-Robinson Oscillatory 
ODE solved by Sunday et al. (2015) 

y' = L(y - sin x) + cos x, y( 0) = 0, L = -1, h = 0.1 

( 21 ) 

Exact solution: y(x) = sin X 


Problem 3: Consider the Growth Model 

A bacteria culture is known to grow at a rate proportional to 
the amount present. After one 

hour, 1000 strands of the bacteria are observed in the 
culture; and after four hours, 3000 strands. 
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Find the number of strands of the bacteria present in the 
culture at time t '. 0 < t < 1 

Let N(t) denote the number of bacteria strands in the culture 
at time t, the initial value 

problem modeling this problem is given by, 

dN 

-= 0.366 N, N( 0) = 694 

dt 

( 22 ) 

The exact solution is given by 

N(t) = 694e 0366f 

Problem 4: Consider the SIR Model 

The SIR model is an epidemiological model that computes 
the theoretical number of people infected with a contagious 
illness in a closed population over time. 

The name of this class of models derives from the fact that 
they involve coupled equations relating the number of 
susceptible people S(t), number of people infected I(t) and 
the number of people who have recovered R(t). This is a 
good and simple model for many infectious diseases 

Including measles, mumps and rubella. It is given by the 
following three coupled equations: 

y(t) = l-0.5e~° 5 ' 

Table 1: The exact solution, computed solution, the absolute error and the time taken for the computation in (seconds) of the new 
method for problem 1 

X-value Exact-solution Computed-solution Error Time/sec 


0.10 

0.9048374180359595 

0.9048374164190154 

1.6169441e-09 

0.205696 

0.20 

0.8187307530779814 

0.8187307469199735 

6.1580079e-09 

0.367746 

0.30 

0.7408182206817172 

0.7408182076888079 

1.2992909e-08 

0.551484 

0.40 

0.6703200460356384 

0.6703200243242226 

2.1711416e-08 

0.702339 

0.50 

0.6065306597126323 

0.6065306278521965 

3.1860436e-08 

0.872062 

0.60 

0.5488116360940253 

0.5488115930362602 

4.3057765e-08 

1.089429 

0.70 

0.4965853037914083 

0.4965852486351736 

5.5156235e-08 

1.311644 

0.80 

0.4493289641172203 

0.4493288965201419 

6.7597078e-08 

1.510724 

0.90 

0.4065696597405978 

0.4065695792037555 

8.0536842e-08 

1.704233 

1.00 

0.3678794411714410 

0.3678793476848027 

9.3486638e-08 

1.928272 


^- = jui\-S)-piS 
dt 

^r=-fii-ri+pis 

dt 

dR _ , 

- = -juR + yI 
dt 


(23) 


where /U,y,f3 are positive parameters to be determined. 
Define y to be, 

y = S + / + R 

and adding the equations in (23) above, we obtain the 
following evolution equation for y, 

y' = ]iil-y) 

Taking fl = 0.5 and attaching an initial condition 
y(0) = 0.5 (for a particular closed population), we obtain, 

-j- = 0.5(1- t), t(0) = 0.5 

dt 


whose Exact solution is 
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Table 2: The exact solution, computed solution, the absolute error and the time taken for the computation in (seconds) of the new 
method for problem 2 


X-value Exact-solution 

Computed-solution 

Error 

Time/sec 

0.10 

0.0998334166468282 

0.0998334154028488 

1.2439794e-09 

0.181433 

0.20 

0.1986693307950614 

0.1986693259603136 

4.8347478e-09 

0.380374 

0.30 

0.2955202066613387 

0.2955201961494993 

1.0511839e-08 

0.564763 

0.40 

0.3894183423086477 

0.3894183242933311 

1.8015317e-08 

0.751144 

0.50 

0.4794255386041985 

0.4794255115178666 

2.7086332e-08 

0.952488 

0.60 

0.5646424733950294 

0.5646424359271566 

3.7467873e-08 

1.129044 

0.70 

0.6442176872376839 

0.6442176383320178 

4.8905666e-08 

1.332819 

0.80 

0.7173560908995147 

0.7173560297503054 

6.1149209e-08 

1.620844 

0.90 

0.7833269096274749 

0.7833268356745610 

7.3952914e-08 

1.904968 

1.00 

0.8414709848078880 

0.8414708977305652 

8.7077323e-08 

2.125000 


Table 3: The exact solution, computed solution, the absolute error and the time taken for the computation in (seconds) of the new 
method for problem 3 


X-value 

Exact-solution 

Computed-solution Error 

Time/sec 

0.010 

696.5446939492859400 

696.5446939492180700 

6.7871042e-011 

0.022520 

0.020 

699.0987185430900600 

699.0987185427908300 

2.9922376e-010 

0.050186 

0.030 

701.6621079941421600 

701.6621079934537900 

6.8837380e-010 

0.070106 

0.040 

704.2348966406204900 

704.2348966393841500 

1.2363444e-009 

0.090180 

0.050 

706.8171189466110000 

706.8171189446453600 

1.9656454e-009 

0.110484 

0.060 

709.4088095025695100 

709.4088094997416600 

2.8278464e-009 

0.133450 

0.070 

712.0100030257847300 

712.0100030218745800 

3.9101451e-009 

0.152786 

0.080 

714.6207343608434700 

714.6207343557549600 

5.0885092e-009 

0.175301 

0.090 

717.2410384800974700 

717.2410384736124300 

6.4850383e-009 

0.210087 

0.100 

719.8709504841319800 

719.8709504760998900 

8.0320888e-009 

0.232457 


Table 4: The exact solution, computed solution, the absolute error and the time taken for the computation in (seconds) of the new 
method for problem 4 

X-value Exact-solution Computed-solution Error Time/sec 


0.010 

0.5024937604036588 

0.5024937604048754 

1.2165824e-012 

0.043527 

0.020 

0.5049750831254160 

0.5049750831324521 

7.0361494e-012 

0.048093 

0.030 

0.5074440301984686 

0.5074440302153604 

1.6891821e-011 

0.053913 

0.040 

0.5099006633466223 

0.5099006633774158 

3.0793479e-011 

0.059570 

0.050 

0.5123450439858337 

0.5123450440363059 

5.0472182e-011 

0.063933 
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0.060 

0.5147772332257459 

0.5147772332973700 

7.1624151e-011 

0.080116 

0.070 

0.5171972918712168 

0.5171972919729365 

1.0171974e-010 

0.085281 

0.080 

0.5196052804238385 

0.5196052805535286 

1.2969015e-010 

0.093241 

0.090 

0.5220012590834500 

0.5220012592496057 

1.6615576e-010 

0.097912 

0.100 

0.5243852877496430 

0.5243852879546123 

2.0496926e-010 

0.104638 


VI. SUMMARY AND CONCLUSION 


This paper examined the derivation, analysis and has 
demonstrated a successful implementation of seventh-eight 
step in form of linear multistep method. The method uses 
block approach for the solution of first order differential 

/ t, .11315 

equations. The results for V I X n , . I, / — U, —, —, —, —, — 

• v n+] > 8 4 8 2 8 


7 

and — were obtained in block forms. The new method 

8 

speed up the computational processes, less burden in the 
implementations and also increase the rate of convergence 
of the solutions. The running times for the problem solved 
are less than four seconds as shown in Table 1-4. The exact 
solution of the method is in good agreement with the 
numerical solution. The method was found to be zero stable, 
consistence, convergent and P-stable as show in figure 1. 
Hence we conclude that seventh-eight method is reliable 
and gives good results with lesser time for computation of 
first order ordinary differential equations. 
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